 ####################################
# 13_brazil_competing_principals.R #
####################################


####################################
#### COMPETING PRINCIPALS MODEL ####
####################################
rm(list=ls())
library (runjags)
library(rjags)
library (coda)
library (reshape2)
library (random)
library (scales)
library(car)

# print (version)
# platform       x86_64-apple-darwin17.0     
# arch           x86_64                      
# os             darwin17.0                  
# system         x86_64, darwin17.0          
# status                                     
# major          4                           
# minor          2.3                         
# year           2023                        
# month          03                          
# day            15                          
# svn rev        83980                       
# language       R                           
# version.string R version 4.2.3 (2023-03-15)
# nickname       Shortstop Beagle            


# Set working directory
savePathRep <- "/all_data/"
#savePathRep <- "~/Documents/Research Projects/Party Effects/JOP_replication_to_submit/all_data/"

setwd (savePathRep)

RandomSeed<- read.csv("brazil_seeds_2MPE.csv")


##Load Initial values for the mean of ideal points and item parameters
load("brazil_initial_values.Rda")


for (i in 1:11){
    for (j in 2:3){
        
        if(i==11) sn<- "mean" else sn<- i
        seednumber<- RandomSeed[i,j]
        set.seed (seednumber)
        if (j==2) { version2save <- "a" } else { version2save <- "b" }
        
        load (paste0("brazil_rollcall_long_format_", sn, ".Rdata"))
        
        
        ##Construct additional variables
        Procedure.PT<- Vote.Features$PT.leader.request.procedure.x * PT.dummy*Pressure$pressure
        Procedure.PSDB<- Vote.Features$PSDB.leader.request.procedure.x * PSDB.dummy*Pressure$pressure
        Procedure.DEM<- Vote.Features$DEM.leader.request.procedure.x * DEM.dummy*Pressure$pressure
        Procedure.PMDB<- Vote.Features$PMDB.leader.request.procedure.x * PMDB.dummy*Pressure$pressure
        Procedure.PP<- Vote.Features$PP.leader.request.procedure.x * PP.dummy*Pressure$pressure
        Procedure.PR<- Vote.Features$PR.leader.request.procedure.x * PR.dummy*Pressure$pressure
        Procedure.leader.pressure<- Procedure.PT+Procedure.PSDB+Procedure.DEM+Procedure.PMDB+Procedure.PP+Procedure.PR
        Procedure.leader<- pmax(Vote.Features$PT.leader.request.procedure.x, Vote.Features$PSDB.leader.request.procedure.x, Vote.Features$DEM.leader.request.procedure.x, Vote.Features$PMDB.leader.request.procedure.x, Vote.Features$PP.leader.request.procedure.x, Vote.Features$PR.leader.request.procedure.x)
        
        pcb<- ifelse(Vote.Features$pcb.request.x==1, 1, 0)
        R28 <- PT.dummy*pcb*Pressure$pressure
        R29 <- PSDB.dummy*pcb*Pressure$pressure
        R30 <- DEM.dummy*pcb*Pressure$pressure
        R31 <- PMDB.dummy*pcb*Pressure$pressure
        R32 <- PR.dummy*pcb*Pressure$pressure
        R33 <- PP.dummy*pcb*Pressure$pressure
        
        potential.fix.bill<- c(25, 41, 42,71, 240, 241)
        
        
        compete.bis.v="model {
        for (i in 1:n.obs) {
            probit(CJR[i]) <- beta[vote[i]]*theta[dep[i]] + eta[vote[i]]*delta[dep[i]] - alpha[vote[i]] +
            lambda[1]*P1[obs[i]] +
            lambda[2]*P2[obs[i]] + lambda[3]*P3[obs[i]] + lambda[4]*P4[obs[i]] +
            lambda[5]*P5[obs[i]] + lambda[6]*P6[obs[i]]
            + lambda[7]*P7[obs[i]]
            +lambda[8]*P8[obs[i]] + lambda[9]*P9[obs[i]] + lambda[10]*P10[obs[i]] +
            lambda[11]*P11[obs[i]]+lambda[12]*P12[obs[i]]
            + lambda[13]*P13[obs[i]]
            +lambda[14]*P14[obs[i]] + lambda[15]*P15[obs[i]] + lambda[16]*P16[obs[i]] +
            lambda[17]*P17[obs[i]] +lambda[18]*P18[obs[i]]
            +lambda[19]*P19[obs[i]]
            +lambda[20]*P20[obs[i]]
            +lambda[21]*P21[obs[i]]
            +lambda[22]*P22[obs[i]]
            +lambda[23]*P23[obs[i]]
            +lambda[24]*P24[obs[i]]
            +lambda[25]*P25[obs[i]]
            +lambda[26]*P26[obs[i]]
            
            
            y[i] ~ dcat(Q[i,1:3]);         # code Nay as 1, NA as 2, Aye as 3
            Q[i,1] <- pow(probVoteDisagree*(1-CJR[i]), p1[i]) * pow(probVoteAgree*(1-CJR[i]), p2[i]) * pow(probVoteUncertain*(1-CJR[i]), (1-p1[i]-p2[i]))
            Q[i,2] <- 1 - Q[i,1] - Q[i,3]
            Q[i,3] <- pow(probVoteAgree*CJR[i], p1[i]) * pow(probVoteDisagree*CJR[i], p2[i]) * pow(probVoteUncertain*CJR[i], (1-p1[i]-p2[i]))
            # Try first with a single probVoteAgree and probVoteDisagree; this means that changes of abstaining
            # are assumed identical across individuals within parties, and across parties
        }
        # we separate legislators for which we do not estimate separate party pressures
        # In reading data, make sure that these are sorted toward the end
        # PRIORS
        probVoteDisagree ~ dbeta(1,1)  # This is a flat prior
        probVoteAgree ~ dbeta(1,1)     # This is a flat prior
        probVoteUncertain ~ dbeta(1,1) # This is a flat prior
        for(j in 1:26) { lambda[j] ~ dnorm(0,1) }
        for(j in 1:n.item) { alpha[j] ~ dnorm(0,1) }
        # Beta (discrimination, dimension 1)
        for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0,1) }
        beta[fix.bill[1]] ~ dnorm(0,1)
        for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0,1) }
        beta[fix.bill[2]]  ~ dnorm(0,1) #
        for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { beta[j]  ~ dnorm(0,1) }
        beta[fix.bill[3]] ~ dnorm(0,1)
        for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { beta[j]  ~ dnorm(0,1) }
        beta[fix.bill[4]] <- -1.3
        for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { beta[j]  ~ dnorm(0,1) }
        beta[fix.bill[5]] <- 1.6
        for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { beta[j]  ~ dnorm(0,1) }
        beta[fix.bill[6]] <- -3
        for(j in (fix.bill[6]+1):n.item) { beta[j]  ~ dnorm(0,1) }
        # Eta (discrimination, dimension 2)
        for(j in 1:(fix.bill[1]-1)) { eta[j]  ~ dnorm(0,1) }
        eta[fix.bill[1]] <- -2
        for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { eta[j]  ~ dnorm(0,1) }
        eta[fix.bill[2]] ~ dnorm(0,1)
        for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { eta[j]  ~ dnorm(0,1) }
        eta[fix.bill[3]] ~ dnorm(0,1)
        for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { eta[j]  ~ dnorm(0,1) }
        eta[fix.bill[4]]  <- 2
        for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { eta[j]  ~ dnorm(0,1) }
        eta[fix.bill[5]] <- -1.2#
        for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { eta[j]  ~ dnorm(0,1) }
        eta[fix.bill[6]] ~ dnorm(0,1) #
        for(j in (fix.bill[6]+1):n.item) { eta[j]  ~ dnorm(0,1) }
        # ideal points
        for(i in 1:n.legs)  { theta[i] ~ dnorm(0,1) }
        for(i in 1:n.legs)  { delta[i] ~ dnorm(0,1) }
        ## Notes
        # n.party.legislators, n.legislators, and n.item are data (counters for the size of vectors)
        # y, party, party.position, pressure, and R are all data:
        # y is an i*j rollcall matrix with 1 (Aye), 0 (Nay), and NA (unobserved)
        # party is a vector of length i with 1=PRI, 2=PAN, 3=PRD
        # party.position is a matrix of size 3*j with the parties' official positions (1 or 0; NA not allowed at this point)
        # pressure is a matrix of size 3*j with continuous scores for each party's pressure on vote j
        # R is an i*j indicator matrix for legislators/votes for which pressure is expected to be effective
    }"
    
    
    position.vec<- position.vec.strict
    
    y1=ifelse (molten.roll.call$rc==-9 | (molten.roll.call$rc==-7 & molten.date.roll.call[,14] != "AJ"), 2, ifelse (molten.roll.call$rc==1, 3, ifelse (molten.roll.call$rc==0, 1, NA)))
    
    compete.data <- dump.format(list(y=y1
    , p1=ifelse(position.vec== 1, 1, 0)
    , p2=ifelse(position.vec==-1, 1, 0)
    , n.obs=nrow(molten.roll.call)
    , obs=c(1:nrow(molten.roll.call))
    , vote=as.numeric(molten.roll.call$vote)
    , dep=molten.party.member$final.leg.no
    , n.legs=max(molten.party.member$final.leg.no)
    , n.item=max(as.numeric(molten.roll.call$vote))
    , P1=R40
    , P2=R41
    , P3=R42
    , P4=R43
    , P5=R44
    , P6=R45
    , P7=R2
    , P8=seniority
    ,P9=R3
    , P10=committee.chair
    , P11=R8
    , P12=district.income
    , P13=R138
    , P14=percentile.rank
    , P15=R124
    , P16=log.district.M
    ,P17=R150
    ,P18=molten.roll.call$year
    ,P19=Procedure.leader.pressure
    ,P20=Procedure.leader
    ,P21=R4
    ,P22=minister
    ,P23=Vote.Features$president.position.x*Pressure$pressure
    ,P24= Vote.Features$president.position.x
    ,P25=pcb*Pressure$pressure
    ,P26=pcb
    , fix.bill=potential.fix.bill))
    
    compete.parameters = c("theta", "delta", "alpha", "beta", "eta", "lambda", "deviance"
    , "probVoteDisagree", "probVoteAgree", "probVoteUncertain")
    
    set.seed (seednumber)
    eta1 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
    beta1 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
    eta1[potential.fix.bill[4:5]] <- c(NA, NA)
    eta1[potential.fix.bill[1]] <- NA
    beta1[potential.fix.bill[4]] <- NA
    beta1[potential.fix.bill[5:6]] <- c(NA, NA)
    compete.inits.1 <- function() {
        dump.format(
        list(
        theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
        delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
        alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
        eta = eta1,
        beta = beta1
        , probVoteAgree = rbeta (1,1,1)
        , probVoteDisagree = rbeta (1,1,1)
        , probVoteUncertain = rbeta (1,1,1)
        , lambda = rnorm (26,0,1)
        ,'.RNG.name'="base::Wichmann-Hill"
        ,'.RNG.seed'= seednumber
        ))
    }
    
    
    
    date()
    set.seed (seednumber+1)
    compete.model <- run.jags(
    model=compete.bis.v,
    monitor=compete.parameters,
    data=compete.data,
    inits= list(compete.inits.1()),
    thin=30, burnin=4000, sample=100,
    check.conv=FALSE, plots=FALSE)
    
    date()
    chainsComp <- mcmc.list(list (compete.model$mcmc[[1]]))
    
    # Uncomment next line to save new version
    # save (chainsComp, file=paste0("brazil_competing_final_result_sample_", sn, version2save, ".Rda"))
}
}

